#include "AlbpAnalysis.h"

AlbpAnalysis::AlbpAnalysis(void) 
: m_iNumPoints(8), m_iRadius(1), m_umPatterns(257), m_SidpMatrix(0), m_NeighborSystem(0), 
  m_iNeighbSysWidth(7), m_iNeighbSysHeight(7), m_iTopN(10), m_iRegionSize(64)
{
}

AlbpAnalysis::AlbpAnalysis(int P, int R) 
: m_iNumPoints(P), m_iRadius(R), m_NeighborSystem(0), m_iNeighbSysWidth(7), 
m_iNeighbSysHeight(7), m_iTopN(10), m_SidpMatrix(0), m_iRegionSize(64)
{ 
}

AlbpAnalysis::~AlbpAnalysis(void) 
{
	if(m_SidpMatrix != 0) 
    {
		free(m_SidpMatrix);
	}
}


void AlbpAnalysis::InitAnalyzer(const cv::Mat* pImgData, int pWidthDim, int pNumChannels, int pTopN) 
{
    m_ImgData = pImgData;
    m_iImgWidth = m_ImgData->cols;
    m_iImgHeight = m_ImgData->rows;
    m_iWidthStep = m_ImgData->step;
    m_iNumChannels = pNumChannels;
    m_iTopN = pTopN;
    m_iRegionSize = pWidthDim;
    m_SidpMatrix = new int[m_iRegionSize * m_iRegionSize];
}

/* method analyzeData
 * Given an array data containing w by h elements, perform texture
 * analysis using the Advanced Local Binary Patterns algorithm.
 * 
 * @param seedLoc  location of the seedling within the data array 
 * @param wsize    size in pixels used as the width and height of the 
                   of the seedling region
 * @return featureCounts  a vector of features generated by the ALBP algorithm
 */
std::vector<double>& AlbpAnalysis::AnalyzeData(SeedNode *pSeedLoc, int pWidthDim) 
{
    // Initialize the x & y start coordinates of the area to be analyzed
    m_iRegionStartX = pSeedLoc->pos.x;
    m_iRegionStartY = pSeedLoc->pos.y;

    // Initialize local row and column variables
	int startRow = m_iRegionStartY + m_iRadius, 
        endRow = startRow + pWidthDim - m_iRadius, 
        startCol = m_iRegionStartX + m_iRadius, 
        endCol = startCol + pWidthDim - m_iRadius,
        row, col, label;

	// Run the LBP algorithm over the (w - m_iRadius) by (h - m_iRadius) cropped image 
	// This is done to avoid pixels that will not have enough neighbors to compute LBP
	for(row = startRow; row < endRow ; row++) 
    {
		for(col = startCol; col < endCol; col++) 
        {
			label = AssignLabel(row, col);
			UpdatePattern(AssignGroup(label), col, row);
		}
	}

    // Compute the SIDP Matrix for the input data
    SetNeighborhoodSystem(NULL, 9, 9);
    ComputeSIDP();

	// Compute the GLAM for the input data
	int *glam = new int[4];
	ComputeGLAM(glam);

    // Create the vector containing the features for this region
    double patternCountSum = 0;
    std::vector<std::pair<int, int>> featureHistogram;
    std::tr1::unordered_map<int, ALBPPattern>::iterator pattern;
    for(pattern = m_umPatterns.begin(); pattern != m_umPatterns.end(); ++pattern)
    {
        patternCountSum += pattern->second.getPatternCount();
        featureHistogram.push_back(std::make_pair(pattern->first, pattern->second.getPatternCount()));
    }

    // Sort them by pattern so that when output to a file, all
    // features in the same column correspond to the same pattern.
    sort(featureHistogram.begin(), featureHistogram.end());

    // Create the normalized feature vector
    std::vector<double> featureCounts;
    for(size_t i = 0; i < featureHistogram.size(); i++)
    {
        featureCounts.push_back(featureHistogram[i].second / patternCountSum);
    }

    // Append GLAM measures to feature vector
    featureCounts.push_back(glam[0]);
    featureCounts.push_back(glam[1]);
    featureCounts.push_back(glam[2]);
    featureCounts.push_back(glam[3]);

    // Return the resulting feature vector
    return featureCounts;
}

/* method AssignLabel
 * Given the center pixel location, assignLabel computes the extended
 * Local Binary Pattern (LBP) measure for that pixel according to the
 * Advanced Local Binary Patterns (ALBP) algorithm.
 *
 * @param row row (y) location of the center pixel
 * @param col column (x) location of the center pixel
 */
int AlbpAnalysis::AssignLabel(int pRow, int pCol) 
{
	int label = 0, point; 
	int pixelValue, centerValue = m_ImgData->ptr<uchar>(pRow)[pCol];
	double offset, x, y;

	// Compute the LBP measure, interpolating pixel values as necessary
	for(point = 0; point < m_iNumPoints; point++) 
    {
		// Compute the location of the current neighborhood pixel
		offset = (2.0 * M_PI * point) / m_iNumPoints;
		x = pCol + m_iRadius * cos(offset);
		y = pRow - m_iRadius * sin(offset);

		// Interpolate pixel value
		pixelValue = InterpolateValue(y, x);

		// Compute the sign of difference and threshold the difference using equations
		// 2.8, 2.9 & 2.10, from chapter 2, to compute the final LBP measure (label).
		// Reference: Pietikainen, M., Hadid, A., Zhao, G., Ahonen, T.;
		//   "Computer Vision Using Local Binary Patters (Computational Imaging and Vision)".
		//   2011. Springer-Verlag London Limited. ISBN-10: 0857297473. 
		//   ISBN-10: 978-0857297471.
		if((pixelValue - centerValue) >= 0)
			label |= (1 << point);
	}

	return label;
}

/* method AssignGroup
 * Given the ALBP label, assignGroup computes the group this label belongs
 * to according to the Advanced Local Binary Patterns (ALBP) algorithm. To
 * find the group of a label, the input integer is continually circularly
 * right-shifted until the minimum value for that input integer is found.
 *
 * @param label the ALBP label that will be grouped
 * @return minimum possible value of label found by circular left shift
 */
int AlbpAnalysis::AssignGroup(int pLabel) 
{
	int group = pLabel, i, tempGroup = pLabel;

    // Compute the group for the input label
	for(i = 1; i < m_iNumPoints; i++) 
    {
		tempGroup = CircularLeftShift(tempGroup, m_iNumPoints);
		if(tempGroup < group)
			group = tempGroup;
	}

	return group;
}

/* method CircularLeftShift
 * As the method name suggests, perform a circular left shift operation on
 * the input integer value (mod numBits).
 * 
 * @param value the integer to be circularly right-shifted
 * @param numBits the number of important bits of the input integer
 * @return the result of performing circular left shift on input integer value
 */
int AlbpAnalysis::CircularLeftShift(int pValue, int pNumBits) 
{
	int valLeftShifted = pValue << 1;
	int mask = 1 << pNumBits;
	return (valLeftShifted | ((valLeftShifted & mask) >> pNumBits)) & ~mask;
}

/* method ComputeGLAM
 * Compute the Gray Level Aura Matrix (GLAM) of the Spatial Distribution Information 
 * of Dominant Patterns (SIDP) matrix. A 1 at location (x, y) in the SIDP matrix means
 * there is a dominant pattern there, while a 0 means no dominant pattern exists there.
 * The GLAM measure m(i, j) will be computed for Set i with respect to Set j. That is,
 * m(i,j) = sum(all occurrences of value j in the neighborhood system of value i).
 * 
 * @param pGlam the resulting GLAM computed from the SIDP matrix
 */
void AlbpAnalysis::ComputeGLAM(int* pGlam) 
{
	int row, col;
	// Go through the SIDP matrix and compute the GLAM measure for Set i wrt Set j
	for(row = 0; row < m_iRegionSize; row++) 
    {
		for(col = 0 ; col < m_iRegionSize; col++) 
        {
			// Compute pGlam for set 0 wrt set 0
			pGlam[0] += ComputeGLAMMeasureAt(0, 0, row, col);
			// Compute pGlam for set 0 wrt set 1
			pGlam[1] += ComputeGLAMMeasureAt(0, 1, row, col);
			// Compute pGlam for set 1 wrt set 1
			pGlam[3] += ComputeGLAMMeasureAt(1, 1, row, col);
		}
	}
	// The occurrence of set 1  in the neighborhood system
	// of set 0 are the same as those for set 0 wrt set 1.
	// So, just copy the value after the computation's done.
	pGlam[2] = pGlam[1];
}


/* method ComputeGLAMMeasureAt
 * Compute the Gray Level Aura Matrix (GLAM) measure of the Spatial Distribution  
 * Information of Dominant Patterns (SIDP) matrix between the two sets representing  
 * gray-level i and gray-level j. The GLAM measure will be computed for Set j with 
 * respect to Set i for the location (pCenterX, pCenterY) within the SIDP matrix.
 *
 * @param i gray-level of the first set
 * @param j gray-level of the second set
 * @param pCenterY row (y) location of neighborhood system center pixel in SIDP matrix
 * @param pCenterX column (x) location of neighborhood system center pixel in SIDP matrix
 * @return the pGlam measure for set i with respect to set j using the neighborhood 
 *   structure ns on the SIDP matrix sidp
 */
int AlbpAnalysis::ComputeGLAMMeasureAt(int i, int j, int pCenterY, int pCenterX) 
{
    // The pGlam measure of set i wrt set j in the neighborhood of (pCenterX, pCenterY)
	int measure = 0;

	// As the measure m(Set_i, Set_j) is the number of occurrences of value j
	// in the neighborhood system of all occurrences of value i, return zero
	// if the center pixel is not value i.
	if(m_SidpMatrix[pCenterY * m_iRegionSize + pCenterX] != i)
		return measure;

	// Compute the convolution of the m_iNeighbSysWidth by m_iNeighbSysHeight m_NeighborSystem neighborhood system / structuring 
	// element with the SIDP matrix m_SidpMatrix at the input location (pCenterX, pCenterY)
	int row, col, nsrow = 0, nscol = 0;
	// For every row in the SIDP matrix
	for(row = pCenterY - (m_iNeighbSysHeight >> 1); row < pCenterY + (m_iNeighbSysHeight >> 1); row++) 
    {
		// Ensure the indices are within the height bounds of the SIDP matrix
		if(0 <= row && row <= m_iRegionSize) 
        {
			// For every column in the current row
			for(col = pCenterX - (m_iNeighbSysWidth >> 1); col < pCenterX + (m_iNeighbSysWidth >> 1); col++) 
            {
				// Ensure the indices are within the width bounds of the SIDP matrix
				if(0 <= col && col <= m_iRegionSize) 
                {
					// The GLAM algorithm does not use the center pixel in its computations
					if(col == pCenterX && row == pCenterY) 
                    {
						// Keep track of what column in the neighborhood system that is currently being processed
						nscol++;
						continue;
					}
					// If the neighborhood system value at this pixel location is 1,
					// and the SIDP matrix value == i, then increment the measure.
					if(m_NeighborSystem[nsrow * m_iNeighbSysWidth + nscol] == 1 && m_SidpMatrix[row * m_iRegionSize + col] == j)
						measure++;
				}

				// Keep track of what column in the neighborhood system that is currently being processed
				nscol++;
			}
		}
		// Keep track of what row in the neighborhood system that is currently being processed
		nsrow++;

		// Reset the column
		nscol = 0;
	}

	return measure;
}

/* method ComputeSIDP
 * Using the spatial and frequency information gathered during the extended
 * Local Binary Pattern computation, this method initializes the Spatial
 * Distribution Information of Dominant Patterns (SIDP) matrix. For the top
 * N frequently occurring patterns, a 1 is placed at each pixel location
 * where that pattern occurs while the rest of the matrix is left 0-filled.
 *
 */
void AlbpAnalysis::ComputeSIDP() 
{
    // Reset the SIDP matrix prior to calculation
    memset(m_SidpMatrix, 0, m_iRegionSize * m_iRegionSize * sizeof(int));

	// As we need to sort patterns based on frequency, we'll build
	// a vector containing (group, frequency) pairs.
	std::vector<std::pair<int,int>> sortedPatterns;
	std::tr1::unordered_map <int, ALBPPattern>::iterator pit;
	for(pit = m_umPatterns.begin(); pit != m_umPatterns.end(); ++pit) 
    {
		sortedPatterns.push_back(std::make_pair(pit->second.getPatternCount(), pit->second.getPatternGroup()));
	}

	// Sort the container
	std::sort(sortedPatterns.begin(), sortedPatterns.end(), std::greater_equal<std::pair<int,int>>( ));

	// For the pTopN patterns, put a 1 in the SIDP matrix where they occur
	int i;
	ALBPPattern pattern;
	for(i = 0; i < m_iTopN; i++) 
    {
		// Get the next most frequent pattern
		int patternGroup = sortedPatterns.at(i).second;
		pattern = m_umPatterns[patternGroup];
		std::list<std::pair<int,int>>::const_iterator locIter;

		// Iterate through all of its locations and set those locations in the SIDP matrix to 1
		for(locIter = pattern.patternBegin(); locIter != pattern.patternEnd(); locIter++) 
        {
			const int col = locIter->first - m_iRegionStartX;
			const int row = locIter->second - m_iRegionStartY;
			m_SidpMatrix[row * m_iRegionSize + col] = 1;
		}
	}
}

/* method UpdatePattern
 * Given the pattern key 'group', look in m_umPatterns to determine if 
 * the pattern was previously added. If not found, add pattern, else
 * update the pattern's occurrance count and locations.
 *
 * @param pGroup the pattern id to update
 * @param pRow the pixel row (y) location where this pattern occurs
 * @param pCol the pixel column (x) location where this pattern occurs
 */
void AlbpAnalysis::UpdatePattern(int pGroup, int pCol, int pRow) 
{
	// Create an iterator which will hold the find result
	std::tr1::unordered_map<int, ALBPPattern>::iterator pattern;
	// Look for the input group within the patterns map
	pattern = m_umPatterns.find(pGroup);
	// If not found, add it to the patterns
	if(pattern == m_umPatterns.end()) 
    {
		m_umPatterns.insert(std::make_pair(pGroup, ALBPPattern(pGroup, pCol, pRow)));
	}
	// Otherwise increment its count
	else 
    {
		pattern->second.incrementCount();
		pattern->second.addLocation(std::make_pair(pCol, pRow));
	}
}

/* method InterpolateValue
 * A bilinear interpolation function to determine data values
 * that are not completely aligned within the data grid.
 *
 * @param data the data array
 * @param row the row (y) location of the value that needs to be interpolated
 * @param col the column (x) location of the value that needs to be interpolated
 * @return the value of the interpolated data point
 */
int AlbpAnalysis::InterpolateValue(double pRow, double pCol) 
{
	int rowTop, rowBottom, colLeft, colRight;
	double topLeftPoint, topRightPoint, bottomLeftPoint, bottomRightPoint;
	double topAverage, bottomAverage, colWeight, rowWeight, finalValue;

	// Compute the indices of the nearest neighboring data points
	rowTop = (int) pRow;
	rowBottom = rowTop + 1;
	colLeft = (int) pCol;
	colRight = colLeft + 1;

	// Ensure they are within valid bounds
	if(rowBottom > (m_iImgHeight - 1)) 
    {
		rowTop--;
		rowBottom--;
	}
	if(colRight > (m_iImgWidth - 1)) 
    {
		colLeft--;
		colRight--;
	}

	// Grab the data points needed for interpolation
    const uchar *topRowPtr = m_ImgData->ptr<uchar>(rowTop),
                *btmRowPtr = m_ImgData->ptr<uchar>(rowBottom);
	topLeftPoint = (double) ((unsigned int) topRowPtr[colLeft]);
	topRightPoint = (double) ((unsigned int) topRowPtr[colRight]);
	bottomLeftPoint = (double) ((unsigned int) btmRowPtr[colLeft]);
	bottomRightPoint = (double) ((unsigned int) btmRowPtr[colRight]);

	// Interpolate the values
	colWeight = abs(colLeft - pCol);
	rowWeight = abs(rowTop - pRow);
	topAverage = (1 - colWeight) * topLeftPoint + colWeight * topRightPoint;
	bottomAverage = (1 - colWeight) * bottomLeftPoint + colWeight * bottomRightPoint;
	finalValue = (1 - rowWeight) * topAverage + rowWeight * bottomAverage;
	
	// Round the final value and return it as an integer
	return ((int) (finalValue + 0.5));
}


void AlbpAnalysis::SetNeighborhoodSystem(int* pNeighbSystem, int pNeighbSysWidth, int pNeighbSysHeight) 
{
	m_NeighborSystem = pNeighbSystem;
	m_iNeighbSysWidth = pNeighbSysWidth;
	m_iNeighbSysHeight = pNeighbSysHeight;
    
    if(m_NeighborSystem == NULL)
    {
        m_NeighborSystem = new int[pNeighbSysWidth * pNeighbSysHeight];
        memset(m_NeighborSystem, 1, sizeof(int) * (pNeighbSysWidth * pNeighbSysHeight));
    }
}
